// standard
#include <vector>
#include <string>
#include <iostream>

// boost
#include <boost/format.hpp>
#include <boost/regex.hpp>
#include <boost/lexical_cast.hpp>

// components
#include <portable_timer.h>
#include <print_error.h>
#include <progress_indicator.h>
#include <matrix.h>
#include <matrix_io.h>
#define DEBUG_LEVEL 6
#include <debug_func.h>

#include "get_arg.h"


using boost::format;
using boost::lexical_cast;
using std::vector;
using std::string;
using std::cout;
using std::cin;
using std::cerr;
using namespace bcpp_matrix;


struct t_index_val
{
	unsigned	index;
	float		val;
	t_index_val(unsigned i, float v) : index(i), val(v){}
	bool operator < (const t_index_val& other) const
	{	return index < other.index;}
};


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	get gene in tax2 which is closest to each gene in tax1

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
void get_closest_from_other_species(t_dense_matrix_f& matrix,
									const vector<unsigned>& tax1_genes,
									const vector<unsigned>& tax2_genes,
									vector<t_index_val>& closest_to_tax1)
{
	debug_func_cerr(5);
	typedef vector<unsigned>::const_iterator c_iter;
	for (c_iter ii = tax1_genes.begin(), ii_end = tax1_genes.end(); ii != ii_end; ++ii)
	{
		double minimum = 10000.0;
		unsigned min_index = static_cast<unsigned>(-1);
		for (c_iter jj = tax2_genes.begin(), jj_end = tax2_genes.end(); jj != jj_end; ++jj)
		{
			if (matrix[*ii][*jj] < minimum)
			{
				min_index = *jj;
				minimum = matrix[*ii][*jj];
			}
		}
		closest_to_tax1.push_back(t_index_val(min_index, minimum));
	}
	assert(closest_to_tax1.size() == tax1_genes.size());
}


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	get gene in tax2 which is closest to each gene in tax1

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
void count_closest_from_same_species(t_dense_matrix_f& matrix,
									const vector<unsigned>& tax_genes,
									const vector<t_index_val>& closest_dist,
									vector<t_index_val>& closest_same_tax,
									vector<unsigned>& cnt_closer)
{
	debug_func_cerr(5);
	assert(closest_dist.size() == tax_genes.size());
	typedef vector<unsigned>::const_iterator c_iter;
	for (unsigned i = 0; i < tax_genes.size(); ++i)
	{
		unsigned ii = tax_genes[i];
		unsigned same_count = 0;
		float other_minimum = closest_dist[i].val;

		double same_minimum = other_minimum;
		unsigned same_min_index = static_cast<unsigned>(-1);

		for (c_iter jj = tax_genes.begin(), jj_end = tax_genes.end(); jj != jj_end; ++jj)
		{
			if (matrix[ii][*jj] < other_minimum)
				++same_count;
			if (matrix[ii][*jj] < same_minimum)
			{
				same_min_index = *jj;
				same_minimum = matrix[ii][*jj];
			}
		}
		cnt_closer.push_back(same_count);
		closest_same_tax.push_back(t_index_val(same_min_index, same_minimum));
	}
	assert(cnt_closer.size() == tax_genes.size());
	assert(closest_same_tax.size() == tax_genes.size());

}
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
int main (int argc, char *argv[])
{
	debug_func_cerr(5);
	// we use cout to print progress so flush per call
	std::ios::sync_with_stdio(false);

	t_program_args args;
	try
	{
	if (!process_command_line_options(argc, argv, args))
		return 1;




	//
	//	read data from input file
	//
	t_matrix_f	data;

	args.VERBOSE() << "\tReading matrix...\n";
	{
		t_progress_indicator dots(args.VERBOSE(), 1);
		dots.use_timer();
		unsigned cnt_duplicate_pairs = 0;
		read_matrix(cin,  data,  args.input_format, cnt_duplicate_pairs, dots);
		if (cnt_duplicate_pairs)
			args.VERBOSE() << "\t" << cnt_duplicate_pairs << " duplicate pairs found.\n";
	}


	//
	// divide genes into one species vs. another
	//
	args.VERBOSE() << "\t" << "Identifying genes from each species\n";
	vector<unsigned> tax1_genes;
	vector<unsigned> tax2_genes;
	const boost::regex  tax1_regex(args.taxon1);
	const boost::regex  tax2_regex(args.taxon2);
	for (unsigned i = 0; i < data.names.size(); ++i)
	{
		if (regex_search(data.names[i], tax1_regex, boost::match_default))
			tax1_genes.push_back(i);
		else if (regex_search(data.names[i], tax2_regex, boost::match_default))
			tax2_genes.push_back(i);
		else
			args.VERBOSE() << "Error:\n\t" << data.names[i] << "matches neither species 1 ("
							<< args.taxon1 << ") or species 2 (" << args.taxon2 << ")\n";
	}

	//
	// closest ortholog comparison
	//
	args.VERBOSE() << "\t" << "Closest orthologues to species 1\n";
	vector<t_index_val> closest_to_tax1;
	get_closest_from_other_species(data.matrix, tax1_genes, tax2_genes, closest_to_tax1);

	args.VERBOSE() << "\t" << "Closest orthologues to species 2\n";
	vector<t_index_val> closest_to_tax2;
	get_closest_from_other_species(data.matrix, tax2_genes, tax1_genes, closest_to_tax2);


	//
	// count species specific duplication
	//
	args.VERBOSE() << "\t" << "Apparent Species 1 specific paralogs\n";
	vector<unsigned> cnt_closer1;
	vector<t_index_val> closest_same_tax1;
	count_closest_from_same_species(data.matrix, tax1_genes, closest_to_tax1,
											closest_same_tax1, cnt_closer1);


	args.VERBOSE() << "\t" << "Apparent Species 2 specific paralogs\n";
	vector<unsigned> cnt_closer2;
	vector<t_index_val> closest_same_tax2;
	count_closest_from_same_species(data.matrix, tax2_genes, closest_to_tax2,
											closest_same_tax2, cnt_closer2);

	//
	// output results for species 1
	//
	if (args.header)
		cout << "name\tclosest_other\tdistance\tcount_same_closer\tclosest_same\tdistance\n"
				"----\t-------------\t--------\t-----------------\t------------\t--------\n";

	for (unsigned i = 0; i < closest_to_tax1.size(); ++i)
	{
		cout		<< data.names[tax1_genes[i]] << "\t";
		if (closest_to_tax1[i].index != static_cast<unsigned>(-1))
			cout	<< data.names[closest_to_tax1[i].index] << "\t"
					<< closest_to_tax1[i].val << "\t";
		else
			cout	<< "\t\t";
		cout		<< cnt_closer1[i] << "\t";
		if (closest_same_tax1[i].index != static_cast<unsigned>(-1))
			cout	<< data.names[closest_same_tax1[i].index] << "\t"
					<< closest_same_tax1[i].val << "\n";
		else
			cout	<< "\t\n";

	}

	for (unsigned i = 0; i < closest_to_tax2.size(); ++i)
	{
		cout		<< data.names[tax2_genes[i]] << "\t";
		if (closest_to_tax2[i].index != static_cast<unsigned>(-1))
			cout	<< data.names[closest_to_tax2[i].index] << "\t"
					<< closest_to_tax2[i].val << "\t";
		else
			cout	<< "\t\t";
		cout		<< cnt_closer2[i] << "\t";
		if (closest_same_tax2[i].index != static_cast<unsigned>(-1))
			cout	<< data.names[closest_same_tax2[i].index] << "\t"
					<< closest_same_tax2[i].val << "\n";
		else
			cout	<< "\t\n";

	}

	args.VERBOSE() << "\n\tAll finished!!\n\n";

	}
	catch(std::exception& e)
	{
		std_print_error(args.CERR(), e.what());
		return false;
	}

	return 0;
}


